Quantifying the Effects of Topography on Forest Structure in a Temperate Coniferous Montane Forest
Alex Gould
Hypotheses
H1. Forest structure is negatively correlated with increasing
elevation. H2. Forest structure is negatively correlated with increasing
slope angle. H3. Forest structure is negatively correlated with north
and east aspects. H4. Forest structure is controlled/determined by
topographic position index/ landform classifications.
(forest structural homogeneity is positively correlated with increasing
spatial resolution of the topographic position index) H5. Forest
structure is negatively correlated with increasing Topographic
Roughness. H6. Forest structure will be positively correlated with land
cover and forest type.
Site Description
The location chosen for this project is Bighorn National Forest in
Wyoming, United States. Bighorn National Forest is approximately 128
meters long and 48 meters wide, encompassing over 404685 hectares of
land (USDA Forest Service 2021). The elevation of Bighorn National
Forest ranges from 1524 meters, to its highest elevation of 4013 meters
at Cloud Peak (USDA Forest Service 2021). Approximately 66% of the land
cover within Bighorn National Forest is forested and made up of
Douglas-fir, ponderosa pine, lodgepole pine, and spruce (Witt 2008). The
sites chosen for this study are located within the National Elevation
Dataset N35 W108, as shown with the black outline in Figure X, which is
approximately 4828 square kilometers. 
Site Selection
From August 12-21, 2021, seventy-four plots at 10x10 meters were
selected across the study area of Bighorn National Forest. The sites
chosen were representative of the adjacent area, with their suitability
primarily derived by their elevation to ensure they were within the
montane biogeographic zone (between 1700 and 3000 meters). [Elevations
lower than 1700 meters, and greater than 3000 meters were derived from a
1/3 arc-second digital elevation model (DEM) were contained in polygons
and excluded from site selection consideration. The secondary
consideration for site selection was given to the accessibility of the
site. The majority of the sites chosen were on average 100 meters from a
secondary road (dirt or gravel and not maintained), with few sites being
greater than 100 meters from a primary road (paved and maintained).
Topographic and geographic features were also included in the analysis
of site accessibility, with steep inclines or declines, cliffs, major
water body crossings and unstable or hazardous features generally
excluded from consideration. Sites that had obvious signs of logging
activity or forest fires were also excluded from consideration. The last
criteria for potential sites were that they be analogous to their
surrounding area, and contained at least one mature tree species
(Diameter Breast Height (DBH) must be 10 centimeters or greater). After
the previous criteria were met, final site selection was derived
stochastically by walking a number of paces in one direction determined
by a random number generator. The location produced from the previous
step was used as the Northwest coordinate of the plot. 
Field Data:
For each site, a 10m x 10m plot was demarcated. The plots were aligned north-south, east-west, with markers placed at the northwest, northeast, southeast and southwest corners. The preliminary elevation and photos of the plot were taken at the northwest corner (Fig. 3). At each corner, as well as at the center of each plot, a CI-110 Plant Canopy Imager recorded the coordinates, as and a wide-angle image of the plot canopy. Additionally, the Plant Canopy Imager generated canopy data for the plot including the Gap Fraction, Photosynthetically Active Radiation (PAR) Leaf Area Index (LAI), PAR Average, Sunflecks, Canopy Density, and Leaf Angle. Geographic and topographic data about each plot was noted as well as the presence or absence of rocks, boulders, downed woody debris (DWD) or standing dead trees (snags), water and ground junipers. The vegetative species composition of the plot was identified and for trees with a DBH of 10cm or greater, DBH, richness and occurrence was documented.
Libraries
library(moments)
library(gvlma)
library(raster)
library(lmtest)
library(car)
library(leaps)
library(tidyverse)
library(leaflet)
library(ggplot2)
library(ggthemes)
library(tidyr)
library(broom)
library(sp)
library(sf)
library(vegan)
library(dplyr)
library(spData)
library(stats)
library(evaluate)
library(spatialEco)
library(knitr)
library(insol)
library(tmap)
library(geodiv)
library(GmAMisc)
library(viridis)
library(rgdal)
library(lidR)
library(lattice)
library(rasterVis)
library(RColorBrewer)
library(ggcorrplot)
library(cowplot)
library(googleway)
library(ggrepel)
library(ggspatial)
library(rnaturalearth)
library(rnaturalearthdata)
library(citation)
library(MASS)
library(lme4)
knitr::opts_chunk$set(cache=TRUE, warning=FALSE, message=FALSE) Read in Field Data
CanopyData <- read.csv("proj_data/CanopyData.csv")
PlotData <- read.csv("proj_data/PlotData.csv")
TreeData <- read.csv("proj_data/TreeData.csv")
TreeHeightsandAge <-read.csv("proj_data/Plot_tree_stats.csv")
DBHsAndDensity <-read.csv("proj_data/DBH_Density.csv") Bighorn National Forest Shape File
nf <- read_sf("proj_data/S_USA.AdministrativeForest/S_USA.AdministrativeForest.shp") %>%
st_transform(.,crs=32613)
bighorn <- dplyr::filter(nf,FORESTNAME=="Bighorn National Forest") Load DEM Data
DEM <-raster("proj_data/DEM.tif") Elevation of Bighorn National Forest
Elevation of Bighorn National Forest.
Load in Wyoming Landuse Landcover Data
#Read in the LULC Raster
lulc=raster("proj_data/lulc.tif")
# LULC Names
lulc_names =read_csv("proj_data/LULC.csv") %>%
dplyr::select(class=ECOLSYS_LU...2, ID=Value)
# Convert to a Raster
lulc=as.factor(lulc)
# Update the RAT with a Left Join
levels(lulc)=left_join(levels(lulc)[[1]],lulc_names)
#table(values(lulc))Plot of Land Use Land Cover
lulc_colors=data.frame(class=levels(lulc)[[1]]$class)
lulc_colors$col=c("#EBA998","#FFFB00", "#4080BF", "#989F60", "#AF5063", "#3200FF", "#A289A8", "#BE20DF", "#971452", "#6250AF" , "#8EB880", "#624074", "#A0B232", "#B232A0", "#FF00E8","#F1197C", "#7C2B4B", "#7C2B72", "#A0D092", "#D095DE", "#70808F", "#10D4EF", "#AF50A7", "#EF103D","#4F30CF", "#FFC300", "#BCFF00", "#00BBFF", "#2080DF", "#DF8E20", "#9A8A86", "#CF3050", "#708F75", "#8F708D", "#95BCDE", "#8F8F70", "#FFCD00", "#EB98E1", "#543E6C", "#DFBA20", "#FF9300", "#00FF43", "#B25932", "#B2327F", "#543B1B", "#1B543B", "#1B3D54", "#805C73", "#B62B85", "#F119B2", "#F18E19", "#591AAF", "#E0FF00", "#E0FF90")
lulcplot<-rasterVis::gplot(lulc)+
geom_raster(aes(fill=as.factor(value)))+
scale_fill_manual(labels=lulc_colors$class,
values=lulc_colors$col,
name="Landcover Type")+
coord_equal()+
theme(legend.position = "bottom", legend.key.size=unit(.5, "line"), legend.text=element_text(size=8))+
guides(fill=guide_legend(ncol=3,byrow=TRUE))
lulcplot
Rename, Summarize, and Group Data
canopy <- CanopyData %>%
group_by(Plot_Number) %>%
dplyr::select(-Sunflecks) %>%
summarise_if(is.numeric, list(mean = function(x) mean(x, na.rm = TRUE),
sd = function(x) sd(x, na.rm = TRUE)))
plots <- PlotData %>%
dplyr::group_by(Plot_Number) %>%
dplyr::summarize(lat=mean(Lat), lon=mean(Lon))
trees<- TreeData %>%
summarize(Plot_Number, Tree_ID, Species_Scientific_Name, DBH, Aprox_Tree_Age) %>%
dplyr:: group_by(Plot_Number) Study Site Plots
#GGPlot of Study Site
studysiteplots<-ggplot(plots)+
geom_point(mapping=aes(x=lon, y=lat, color=Plot_Number)) +
scale_color_viridis_c(option="turbo")+
labs(x = "Longitude", y ="Latitude",
title ="Study Site Plots", color ="Plot Number") +
theme(plot.title = element_text(hjust = 0.5))
studysiteplots
#Isolate standard lidar study plots
Lidarplots <-plots[-c(7,12,13, 23,24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 43,
51,57, 58,59, 65,66, 69),]
#GGPlot of Lidar Plots
lidarplotsgg<-ggplot(Lidarplots)+
geom_point(mapping=aes(x=lon, y=lat, color=Plot_Number))+
scale_color_viridis_c(option="turbo")+
labs(x = "Longitude", y ="Latitude",
title ="Lidar Plots", color ="Plot Number") +
theme(plot.title = element_text(hjust = 0.5))
lidarplotsgg
Table of Species by Plot
tree1<-trees %>%
dplyr::select(-Tree_ID, -DBH) %>%
dplyr::group_by(Plot_Number,Species_Scientific_Name) %>%
dplyr::summarise(count=n()) %>%
tidyr::spread(Species_Scientific_Name,count) %>%
dplyr::mutate_all(function(x) ifelse(is.na(x),0,x))
treecounts <-table(trees$Species_Scientific_Name)
view(tree1)
view(treecounts)
tree_coords <- merge(trees, plots)Trees
This graph displays the location of species across plots.
spatial_tc<-ggplot(data=tree_coords, mapping=aes(x=lon, y=lat, color=Species_Scientific_Name))+
geom_point()+
scale_color_viridis_d(option="turbo")+
labs(x = "Latitude", y = "Longitude", title ="Spatial Occurrence of Tree Species", color ="Species Scientific Name") +
theme(plot.title = element_text(hjust = 0.5))
spatial_tc
Species Occurence Throughout Plots
This bar graph displays the total observations of each species across all plots. The Pinus contorta or Lodgepole Pine was observed the most frequently with more than 350 sightings. The Pinus flexilis, Limber Pine, and the Pinus ponderosa, the Ponderosa Pine each only had one sighting across the plots.
tree_counts_bargraph<-barplot(treecounts,
main= "Species Occurence Throughout Plots",
horiz=TRUE,
xlab="Occurence",
ylab= "Species",
col = c("#DD8317", "#3A9B44", "#47ACC2", "#EACF4F", "#EA5F4f", "#DD1794","#FFC300"),
names.arg=c("Juniperus occidentalis", "Picea engelmannii","Pinus contorta", "Pinus ponderosa", "Pinus flexilis", "Populus tremuloides", "Pseudotsuga menziesii"),
cex.names=0.2, legend=TRUE)
Calculate TPI, Slope angle, Slope aspect and Elevation
tpi10m <- raster("proj_data/tpi10m.tif")
tpi100m <- raster("proj_data/tpi100m.tif")
tpi1000m <- raster("proj_data/tpi1000m.tif")
TRI <- raster("proj_data/TRI.tif")
slope <- raster("proj_data/slope.tif")
Slope Angle of Bighorn National Forest.
aspect <- raster("proj_data/aspect.tif")
Slope Aspect of Bighorn National Forest.
eastwest <- raster("proj_data/eastwest.tif")
northsouth <- raster("proj_data/northsouth.tif")Create Landforms at 10 Meter Resolution
SD10m <- sd(tpi10m[],na.rm=T)
SD10m## [1] 0.4458888landform_sd10m <- reclassify(tpi10m, matrix(c(-Inf, -SD10m, 1,
-SD10m, -SD10m/2, 2,
-SD10m/2, 0, 3,
0, SD10m/2, 4,
SD10m/2, SD10m, 5,
SD10m, Inf, 6),
ncol = 3, byrow = T),
right = T)
# Turn it into categorical raster
landform_sd10m <- as.factor(landform_sd10m)
rat_sd10m <- levels(landform_sd10m)[[1]]
rat_sd10m[["landform_sd10m"]] <- c('Valley', 'Lower Slope',
'Flat Area','Middle Slope',
'Upper Slope', 'Ridge')
levels(landform_sd10m) <- rat_sd10m
# Plot the classification
tpi10mplot <-rasterVis::levelplot(landform_sd10m, col.regions = rev(brewer.pal(6,'RdYlBu')),
labels = rat_sd10m$landcover,
main = "10M TPI Landform Classification",
colorkey=list(labels=list(at=1:6, labels=rat_sd10m[["landform_sd10m"]])))
tpi10mplot
histtpi10m <-hist(tpi10m[])
Create Landforms at 100 Meter Resolution
SD100m <- sd(tpi100m[],na.rm=T)
SD100m ## [1] 2.491715# Make Landform Classifications
#Morphologic class De Reu et al. 2013; Weiss (2001)
landform_sd100m <- reclassify(tpi100m, matrix(c(-Inf, -SD100m, 1,
-SD100m, -SD100m/2, 2,
-SD100m/2, 0, 3,
0, SD100m/2, 4,
SD100m/2, SD100m, 5,
SD100m, Inf, 6),
ncol = 3, byrow = T),
right = T)
# Turn it into categorical raster
landform_sd100m <- as.factor(landform_sd100m)
rat_sd100m <- levels(landform_sd100m)[[1]]
rat_sd100m[["landform_sd100m"]] <- c('Valley', 'Lower Slope',
'Flat Area','Middle Slope',
'Upper Slope', 'Ridge')
levels(landform_sd100m) <- rat_sd100m
#writeRaster(landform_sd100m, file="proj_data/landform_sd100m.grd", overwrite=TRUE)
# Plot the classification
tpi100mplot<- rasterVis::levelplot(landform_sd100m, col.regions = rev(brewer.pal(6,'RdYlBu')),
labels = rat_sd100m$landcover,
main = "100m TPI Landform Classification",
colorkey=list(labels=list(at=1:6, labels=rat_sd100m[["landform_sd100m"]])))
tpi100mplot
histtpi100m <- histtpi100<-hist(tpi100m[])
Create Landforms at 1000 Meter Resolution
SD1000m <- sd(tpi1000m[],na.rm=T)
SD1000m ## [1] 29.7036# Make landform classes
#Morphologic class De Reu et al. 2013; Weiss (2001)
landform_sd1000m <- reclassify(tpi1000m, matrix(c(-Inf, -SD1000m, 1,
-SD1000m, -SD1000m/2, 2,
-SD1000m/2, 0, 3,
0, SD1000m/2, 4,
SD1000m/2, SD1000m, 5,
SD1000m, Inf, 6),
ncol = 3, byrow = T),
right = T)
# Turn it into categorical raster
landform_sd1000m <- as.factor(landform_sd1000m)
rat_sd1000m <- levels(landform_sd1000m)[[1]]
rat_sd1000m[["landform_sd1000m"]] <- c('Valley', 'Lower Slope',
'Flat Area','Middle Slope',
'Upper Slope', 'Ridge')
levels(landform_sd1000m) <- rat_sd1000m
#writeRaster(landform_sd1000m, file="proj_data/landform_sd1000m.grd", overwrite= TRUE) #tif wont work for as.factor
# Plot the classification
tpi1000mplot<- levelplot(landform_sd1000m, col.regions = rev(brewer.pal(6,'RdYlBu')),
labels = rat_sd1000m$landcover,
main = "1000m TPI Landform Classification",
colorkey=list(labels=list(at=1:6, labels=rat_sd1000m[["landform_sd1000m"]])))
tpi1000mplot
histtpi1000m<-hist(tpi1000m[])
PlotLandformClassification <-read.csv("proj_data/TPI Counts.csv")Lidar Images of Plots
#### Read in Plots
plotsf<- plots %>%
st_as_sf(coords = c(x="lon", y="lat")) %>%
st_set_crs(.,value=4326)
# Obtain Plot Coordinates in Lidar Projection
lasfile="proj_data/LAS/plot_52_USGS_LPC_WY_Sheridan_2020_D20_13TCK1344.laz"
lproj=readLASheader(lasfile) %>%
crs()
#10m polys
plot_polys <- plotsf %>%
st_transform(.,crs="+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m no_defs") %>%
st_buffer(dist = 50,endCapStyle="SQUARE") %>%
st_transform(.,crs=st_crs(4326))
# Plot Specific Change
plotgeo <-plot_polys[52,] %>%
st_transform(.,crs=lproj)
# Read in New LAS File
lasPlots <-readLAS(lasfile,
filter = paste("-keep_xy ",paste(st_bbox(plotgeo),
collapse=" "))) %>%
clip_roi(geometry= plotgeo)
plot(lasPlots) 
Calculating Tree Heights
#dtm, dsm, chm
dsm2 <- lidR::grid_canopy(lasPlots, res=2, p2r())
dtm2 <- lidR::grid_terrain(lasPlots, res= 2, algorithm=tin())
#Crop!
chm2<-dsm2-dtm2
chmcrop <- crop(chm2, plotgeo)
plot(chmcrop)
#TRUE CROPPED MEAN
chmmean<- cellStats(chmcrop, stat=mean)
view(chmmean)
#TRUE CROPPED SD
chmstd<- cellStats(chmcrop, stat=sd)
view(chmstd)
#Max Value
chmmax <- cellStats(chmcrop, stat=max)
view(chmmax)
#Min Value
chmmin <- cellStats(chmcrop, stat=min)
view(chmmin)
#Plot and Hist
plot(chmcrop, col = height.colors(25))
hist(chmcrop)
hist(TreeHeightsandAge$Max_Tree_Height)
hist(TreeHeightsandAge$Mean_Tree_Height)
Calculate the Diversity
diversity=data.frame(Plot_Number=tree1$Plot_Number,div=diversity(tree1))Species Diversity
This bar graph shows the diversity within each plot. Plot 2 high the highest diversity at 1.2366849, and plot 60 had the lowest diversity of 0.0836497.
spatial_div<-ggplot(data= diversity, mapping=aes(x=Plot_Number, y=div, fill=div))+
geom_bar(stat = "identity")+
scale_fill_viridis_c(option="turbo")+
labs(x = "Plot Number", y = "Diversity", title ="Species Diversity Across Plots")
spatial_div
Create a Merged Dataset
Merged <- merge(canopy, plots) %>%
merge(DBHsAndDensity)%>%
merge(diversity) %>%
st_as_sf(coords = c(x="lon", y="lat")) %>%
st_set_crs(.,value=4326) %>%
st_transform(.,crs=st_crs(DEM))
Merged_With_TreeHeightsandAge <- merge(canopy, plots) %>%
merge(TreeHeightsandAge) %>%
st_as_sf(coords = c(x="lon", y="lat")) %>%
st_set_crs(.,value=4326) %>%
st_transform(.,crs=st_crs(DEM))Plotted Species Diversity
Plotted_div<-ggplot(data=Merged, mapping=aes(x=plots$lon, y=plots$lat, color =div))+
geom_point()+
scale_color_viridis_c(option="turbo")+
labs(x = "Longitude", y = "Latitude", title ="Species Occurene Throughout Plots", color ="Diversity")
Plotted_div
Stack All Topographic Variables
all_topo <- raster::stack(DEM,
TRI,
slope, aspect,
eastwest, northsouth,
tpi10m, tpi100m, tpi1000m)
plot(all_topo)
Combine and Clean all Data at 10 Meter Resolution
envi <- raster::extract(all_topo, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean, file="proj_data/envi.tif", overwrite=TRUE
) %>%
scale() #this makes the regression coefficients comparable
envi<-as.data.frame(envi)
bind<- bind_cols(Merged,envi)
#LULC Introduction
envi2 <- raster::extract(lulc, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"), small=TRUE, fun=mean, file="proj_data/envi2.tif", overwrite=TRUE
) %>%
as.data.frame() %>%
left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)
bighornData <- bind_cols(bind,as.data.frame(envi2)) %>%
mutate(lulc=as.factor(lulc)) %>%
mutate(lulc_reduced= as.factor(case_when(
lulc==149~149,
lulc==151~151,
lulc==155~155,
TRUE~1,
)))
# dplyr::select(-"X")Statistics Per Plot
This image displays the data collected at each plot in a spatial format. The graphic in the top left shows the plots, numbered. The next image shows the mean Gap Fraction per plot. Followed by the mean PAR LAI per plot, the mean PAR Average per plot and the mean Canopy Density per plot. The second row of images shows the mean Leaf Angle per plot, the standard deviation of the Gap Fraction per plot, the standard deviation of the PAR LAI per plot, the standard deviation of the PAR Average per plot and the standard deviation of the Canopy Density per plot.
plot(bind)
title(main = "Extracting Points", line= 8)
Not Scaled- For Individual Maps!
enviNS <- raster::extract(all_topo, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean, file="proj_data/enviNS.tif", overwrite=TRUE)
enviNS<-as.data.frame(enviNS)
bindNS<- bind_cols(Merged,enviNS)
treesandbindNS<- merge(trees, bindNS)
#Add LULC
envi2NS <- raster::extract(lulc, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"), small=TRUE, fun=mean, file="proj_data/envi2NS.tif", overwrite=TRUE
) %>%
as.data.frame() %>%
left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)
bighornDataNS <- bind_cols(bindNS,as.data.frame(envi2NS)) %>%
mutate(lulc=as.factor(lulc)) %>%
mutate(lulc_reduced= as.factor(case_when(
lulc==149~149,
lulc==151~151,
lulc==155~155,
TRUE~1,
)))Tree Height Combine and Clean at 10 Meter Resolution
enviTree <- raster::extract(all_topo, st_buffer(Merged_With_TreeHeightsandAge,
dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean,
file="proj_data/enviTree.tif", overwrite=TRUE
) %>%
scale() #this makes the regression coefficients comparable
enviTree<-as.data.frame(enviTree)
bindTree<- dplyr::bind_cols(Merged_With_TreeHeightsandAge,enviTree) %>%
dplyr::rename(Average_Stand_Age=Average_Stamd_Age)
# LULC
envi2Tree <- raster::extract(lulc, st_buffer(Merged_With_TreeHeightsandAge, dist = 5, endCapStyle="SQUARE"),
small=TRUE, fun=mean,file="proj_data/envi2Tree.tif", overwrite=TRUE) %>%
as.data.frame() %>%
left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)
bighornDataTreeHeightsAndAge <- bind_cols(bindTree,as.data.frame(envi2Tree)) %>%
mutate(lulc=as.factor(lulc)) %>%
dplyr::select(-"X") %>%
mutate(lulc_reduced= as.factor(case_when(
lulc==149~149,
lulc==151~151,
lulc==155~155,
TRUE~1,
)))Tree Height Not Scaled for Graphing
enviTreeNS <- raster::extract(all_topo, st_buffer(Merged_With_TreeHeightsandAge,
dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean,
file="proj_data/enviTreeNS.tif", overwrite=TRUE
)
enviTreeNS<-as.data.frame(enviTreeNS)
bindTreeNS<- dplyr::bind_cols(Merged_With_TreeHeightsandAge,enviTreeNS) %>%
dplyr::rename(Average_Stand_Age=Average_Stamd_Age)
#LULC
envi2TreeNS <- raster::extract(lulc, st_buffer(Merged_With_TreeHeightsandAge, dist = 5, endCapStyle="SQUARE"), small=TRUE, fun=mean,
file="proj_data/envi2TreeNS.tif", overwrite=TRUE) %>%
as.data.frame() %>%
left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)
NSbighornDataTreeHeightsAndAge <- bind_cols(bindTreeNS,as.data.frame(envi2TreeNS)) %>%
mutate(lulc=as.factor(lulc)) %>%
dplyr::select(-"X") %>%
mutate(lulc_reduced= as.factor(case_when(
lulc==149~149,
lulc==151~151,
lulc==155~155,
TRUE~1,
)))Interactive Study Site Plots
# Interactive Map of All Plots
tmap::tmap_mode(mode='view')
tm_shape(bindTree)+
tm_dots(col='black')#Isolate Interactive Lidar Plots
LidarplotsInt <-bindTree[-c(7,12,13, 23,24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 43, 51,57, 58,59, 65,66, 69),]
#Interactive Map of Lidar Plots
tmap::tmap_mode(mode='view')
tm_shape(LidarplotsInt, )+
tm_dots(col='blue')Forest Structure Correlations at 10 Meters
view(bighornData)
cor_bighorn10m <- bighornData %>% data.frame() %>%
dplyr::select(
#"Plot_Number" ,
"Gap_Fraction_mean",
"Stand_Density",
"Average_DBH",
"DBH_Sd",
"PAR_LAI_mean" ,
"PAR_Average_mean" ,
"Canopy_Density_mean",
"Leaf_Angle_mean" ,
"Gap_Fraction_sd",
"PAR_LAI_sd" ,
"PAR_Average_sd" ,
"Canopy_Density_sd",
"Leaf_Angle_sd" ,
"div",
"DEM" ,
"TRI" ,
"slope" ,
"aspect" ,
"eastwest" ,
"northsouth" ,
"tpi10m" ,
"tpi100m" ,
"tpi1000m" ,
) %>%
na.omit()
my_cor10<-cor(cor_bighorn10m)
view(my_cor10)
my_pmat10<-cor_pmat(cor_bighorn10m)
view(my_pmat10%>% format(scientific=F))
ggcorrplot(corr = my_cor10, p.mat = my_pmat10)
my_cor10kendall<-cor(cor_bighorn10m, method="kendall")
view(my_cor10kendall)
my_pmat10kendall<-cor_pmat(cor_bighorn10m, method="kendall")
view(my_pmat10kendall %>% format(scientific=F))
ggcorrplot(corr = my_cor10kendall, p.mat = my_pmat10kendall)
Tree Height Correlations at 10 Meeters
cor_bighornTreeHeightsAndAge <-bighornDataTreeHeightsAndAge %>%
data.frame() %>%
dplyr::select(
"Mean_Tree_Height" ,
"Sd_Tree_Height",
"Max_Tree_Height",
"Min_Tree_Height",
"Average_Stand_Age",
"DEM" ,
"TRI" ,
"slope" ,
"aspect" ,
"eastwest" ,
"northsouth" ,
"tpi10m" ,
"tpi100m" ,
"tpi1000m" ,
) %>%
na.omit()
my_corTrees<-cor(cor_bighornTreeHeightsAndAge)
my_pmatTrees<-cor_pmat(cor_bighornTreeHeightsAndAge)
ggcorrplot(corr = my_corTrees,p.mat = my_pmatTrees)
my_corTreesKendall<-cor(cor_bighornTreeHeightsAndAge, method="kendall")
my_pmatTreesKendall<-cor_pmat(cor_bighornTreeHeightsAndAge, method="kendall")
ggcorrplot(corr = my_corTreesKendall,p.mat = my_pmatTreesKendall)
Model Building
####Gap Fraction mean ####
lm_GapFractionmMean10_exhaustive<-regsubsets(Gap_Fraction_mean~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+tpi1000m+lulc_reduced, data=bighornData,intercept=TRUE,method="exhaustive", nvmax =12,nbest = 3,)
lm_GapFractionmMean10<-lm_GapFractionmMean10_exhaustive %>%
broom::tidy()
view(lm_GapFractionmMean10)
####Gap Fraction sd ####
lm_GapFraction_sd10m_exhaustive<-regsubsets(Gap_Fraction_sd~DEM+TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+tpi1000m+ lulc_reduced,data=bighornData, intercept=TRUE, method="exhaustive", nvmax =12, nbest = 3)
lm_GapFraction_sd10m<-lm_GapFraction_sd10m_exhaustive %>%
broom::tidy()
view(lm_GapFraction_sd10m)
#####################################################################################
####PAR LAI mean ####
lm_PAR_LAI_mean_10m_exhaustive<-regsubsets(PAR_LAI_mean~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3,)
lm_PAR_LAI_mean_10m<-lm_PAR_LAI_mean_10m_exhaustive %>%
broom::tidy()
view(lm_PAR_LAI_mean_10m)
####PAR LAI sd ####
lm_PAR_LAI_sd_10m_exhaustive<-regsubsets(PAR_LAI_sd~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_PAR_LAI_sd_10m<-lm_PAR_LAI_sd_10m_exhaustive %>%
broom::tidy()
view(lm_PAR_LAI_sd_10m)
#####################################################################################
####PAR Average mean ####
lm_PAR_Average_mean_exhaustive<-regsubsets(PAR_Average_mean~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_PAR_Average_mean_10m<-lm_PAR_Average_mean_exhaustive %>%
broom::tidy()
view(lm_PAR_Average_mean_10m)
####PAR Average SD ####
lm_PAR_Average_sd10m_exhaustive<-regsubsets(PAR_Average_sd~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_PAR_Average_sd10m<-lm_PAR_Average_sd10m_exhaustive %>%
broom::tidy()
view(lm_PAR_Average_sd10m)
#####################################################################################
####Canopy Density Mean ####
lm_Canopy_Density_mean_10m_exhaustive <-regsubsets(Canopy_Density_mean~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_Canopy_Density_mean_10m<-lm_Canopy_Density_mean_10m_exhaustive %>%
broom::tidy()
view(lm_Canopy_Density_mean_10m)
####Canopy Density SD####
lm_Canopy_Density_sd_10m_exhaustive <-regsubsets(Canopy_Density_sd~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_Canopy_Density_sd_10m<-lm_Canopy_Density_sd_10m_exhaustive %>%
broom::tidy()
view(lm_Canopy_Density_sd_10m)
#####################################################################################
#### Leaf Angle mean####
lm_Leaf_Angle_mean_10_exhaustive <-regsubsets(Leaf_Angle_mean~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_Leaf_Angle_mean_10<-lm_Leaf_Angle_mean_10_exhaustive %>%
broom::tidy()
view(lm_Leaf_Angle_mean_10)
####Leaf Angle sd ####
lm_Leaf_Angle_sd_10_exhaustive <-regsubsets(Leaf_Angle_sd~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_Leaf_Angle_sd_10<-lm_Leaf_Angle_sd_10_exhaustive %>%
broom::tidy()
view(lm_Leaf_Angle_sd_10)
#####################################################################################
####Diversity ####
lm_Diversity_exhaustive <-regsubsets(div~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_Diversity_10m<-lm_Diversity_exhaustive %>%
broom::tidy()
view(lm_Diversity_10m)
#####################################################################################
####Average DBH ####
lm_DBH_mean_10m_exhaustive <-regsubsets(Average_DBH~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_DBH_mean_10m<-lm_DBH_mean_10m_exhaustive %>%
broom::tidy()
view(lm_DBH_mean_10m)
#####################################################################################
####Dbh SD ####
lm_DBH_sd_10m_exhaustive <-regsubsets(DBH_Sd~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_DBH_sd_10m<-lm_DBH_sd_10m_exhaustive %>%
broom::tidy()
view(lm_DBH_sd_10m)
#####################################################################################
####Stand Density ####
lm_Stand_Density_10m_exhaustive <-regsubsets(Stand_Density~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest = 3)
lm_Stand_Density_10m<-lm_Stand_Density_10m_exhaustive %>%
broom::tidy()
view(lm_Stand_Density_10m)
###################################################################################
####Mean Tree Height ####
lm_TreeHeightMean_exhaustive<-regsubsets(Mean_Tree_Height~DEM+TRI+slope+aspect+ eastwest+northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced, data=bighornDataTreeHeightsAndAge, intercept=TRUE, method="exhaustive", nvmax =10, nbest = 3)
lm_TreeHeightMean<-lm_TreeHeightMean_exhaustive %>%
broom::tidy()
view(lm_TreeHeightMean)
###################################################################################
####SD Tree Height ####
lm_TreeHeight_SD_exhaustive<-regsubsets(Sd_Tree_Height~DEM+TRI+slope+aspect+ eastwest+northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced, data=bighornDataTreeHeightsAndAge, intercept=TRUE, method="exhaustive", nvmax =10, nbest = 3)
lm_TreeHeightSD<-lm_TreeHeight_SD_exhaustive %>%
broom::tidy()
view(lm_TreeHeightSD)
##################################################################################
#Max Tree Height
lm_Max_TreeHeight_exhaustive<-regsubsets(Max_Tree_Height~DEM+TRI+slope+aspect+ eastwest+northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced, data=bighornDataTreeHeightsAndAge, intercept=TRUE, method="exhaustive", nvmax =10, nbest = 3)
lm_Max_TreeHeight<-lm_Max_TreeHeight_exhaustive %>%
broom::tidy()
view(lm_Max_TreeHeight)
##################################################################################
#Average Stand Age (tpi100 and 1000 individually, and Landforms100 and 1000)
lm_Average_Stand_Age10m_exhaustive<-regsubsets(Average_Stand_Age~DEM+TRI+slope+aspect+ eastwest+northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced, data=bighornDataTreeHeightsAndAge, intercept=TRUE, method="exhaustive", nvmax =10, nbest = 3)
lm_Average_Stand_Age10m<-lm_Average_Stand_Age10m_exhaustive %>%
broom::tidy()
view(lm_Average_Stand_Age10m)Forest Structure Regressions at 10 Meters
#Gap Fraction mean
#p-value: 2.378e-05
GapFraction_mean_10m<-lm(Gap_Fraction_mean~ TRI+ tpi1000m+lulc_reduced,
data=bighornData)
summary(GapFraction_mean_10m)##
## Call:
## lm(formula = Gap_Fraction_mean ~ TRI + tpi1000m + lulc_reduced,
## data = bighornData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.90669 -0.16654 0.00984 0.20381 0.72175
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.11244 0.09494 11.717 < 2e-16 ***
## TRI 0.07958 0.03974 2.003 0.0493 *
## tpi1000m 0.08333 0.04003 2.082 0.0412 *
## lulc_reduced149 0.49340 0.10851 4.547 2.39e-05 ***
## lulc_reduced151 0.60939 0.14439 4.220 7.61e-05 ***
## lulc_reduced155 0.16545 0.15041 1.100 0.2753
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3257 on 66 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.349, Adjusted R-squared: 0.2997
## F-statistic: 7.076 on 5 and 66 DF, p-value: 2.378e-05broom::tidy(GapFraction_mean_10m)## # A tibble: 6 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1.11 0.0949 11.7 8.91e-18
## 2 TRI 0.0796 0.0397 2.00 4.93e- 2
## 3 tpi1000m 0.0833 0.0400 2.08 4.12e- 2
## 4 lulc_reduced149 0.493 0.109 4.55 2.39e- 5
## 5 lulc_reduced151 0.609 0.144 4.22 7.61e- 5
## 6 lulc_reduced155 0.165 0.150 1.10 2.75e- 1broom::glance(GapFraction_mean_10m)## # A tibble: 1 × 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.349 0.300 0.326 7.08 0.0000238 5 -18.3 50.5 66.5
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>broom::augment(GapFraction_mean_10m)## # A tibble: 72 × 11
## .rownames Gap_Fraction_mean TRI tpi1000m lulc_reduced .fitted .resid
## <chr> <dbl> <dbl> <dbl> <fct> <dbl> <dbl>
## 1 1 1.22 -0.973 0.920 1 1.11 0.108
## 2 2 1.72 -1.41 0.545 149 1.54 0.181
## 3 3 1.72 -0.588 0.771 149 1.62 0.0967
## 4 4 1.7 -1.08 1.78 149 1.67 0.0315
## 5 5 1.78 2.53 1.24 149 1.91 -0.130
## 6 6 1.78 -0.944 0.349 149 1.56 0.220
## 7 8 1.72 -1.07 -0.0324 149 1.52 0.202
## 8 9 1.56 -0.182 -0.0256 149 1.59 -0.0292
## 9 10 1.46 0.433 0.687 149 1.70 -0.238
## 10 11 1.58 0.0482 0.117 149 1.62 -0.0394
## # … with 62 more rows, and 4 more variables: .hat <dbl>, .sigma <dbl>,
## # .cooksd <dbl>, .std.resid <dbl>densityPlot(GapFraction_mean_10m$residuals)
AIC(GapFraction_mean_10m)## [1] 50.52386BIC(GapFraction_mean_10m)## [1] 66.46053#
# MVGapFraction_mean_10m<-glm(Gap_Fraction_mean~TRI+ tpi1000m+lulc_reduced,
# family=gaussian(link="log"),
# data=bighornData)
# summary(MVGapFraction_mean_10m)
#BIC(GapFraction_mean_10m, MVGapFraction_mean_10m)
#Gap Fraction sd
#p-value: 0.05405
# GapFraction_sd_10m<-glm(Gap_Fraction_sd~ DEM,
# family = gaussian(),
# data=bighornData)
# summary(GapFraction_sd_10m)
#densityPlot(GapFraction_sd_10m$residuals)
MVGapFraction_sd_10m<-glm(Gap_Fraction_sd~DEM,
family=gaussian(link="log"),
data=bighornData)
summary(MVGapFraction_sd_10m)##
## Call:
## glm(formula = Gap_Fraction_sd ~ DEM, family = gaussian(link = "log"),
## data = bighornData)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.22114 -0.13160 -0.06747 0.07563 0.78255
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.24112 0.07891 -15.728 <2e-16 ***
## DEM -0.14521 0.06412 -2.265 0.0266 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.03671833)
##
## Null deviance: 2.7575 on 72 degrees of freedom
## Residual deviance: 2.6070 on 71 degrees of freedom
## AIC: -30.09
##
## Number of Fisher Scoring iterations: 6AIC(MVGapFraction_sd_10m)## [1] -30.09BIC(MVGapFraction_sd_10m)## [1] -23.21862broom::tidy(MVGapFraction_sd_10m) ## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -1.24 0.0789 -15.7 7.82e-25
## 2 DEM -0.145 0.0641 -2.26 2.66e- 2broom::glance(MVGapFraction_sd_10m)## # A tibble: 1 × 8
## null.deviance df.null logLik AIC BIC deviance df.residual nobs
## <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
## 1 2.76 72 18.0 -30.1 -23.2 2.61 71 73densityPlot(MVGapFraction_sd_10m$residuals)
#####################################################################################
####PAR LAI mean ####
# #p-value: 0.1319
# PAR_LAI_mean_10m<-lm(PAR_LAI_mean~ DEM+tpi10m+lulc_reduced,
# data=bighornData)
# summary(PAR_LAI_mean_10m)
mvPAR_LAI_mean_10m<-glm(PAR_LAI_mean~ DEM+tpi10m+lulc_reduced,
family=gaussian(link="log"),
data=bighornData)
summary(mvPAR_LAI_mean_10m)##
## Call:
## glm(formula = PAR_LAI_mean ~ DEM + tpi10m + lulc_reduced, family = gaussian(link = "log"),
## data = bighornData)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4816 -0.9334 -0.2222 0.9833 4.3572
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.72732 0.20962 3.470 0.000916 ***
## DEM 0.17283 0.08178 2.113 0.038297 *
## tpi10m 0.12018 0.08230 1.460 0.148864
## lulc_reduced149 0.22686 0.22516 1.008 0.317289
## lulc_reduced151 0.16763 0.27760 0.604 0.547980
## lulc_reduced155 0.60219 0.26272 2.292 0.025045 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 2.234312)
##
## Null deviance: 172.23 on 72 degrees of freedom
## Residual deviance: 149.70 on 67 degrees of freedom
## AIC: 273.59
##
## Number of Fisher Scoring iterations: 8densityPlot(mvPAR_LAI_mean_10m$residuals)
AIC(mvPAR_LAI_mean_10m)## [1] 273.5912BIC(mvPAR_LAI_mean_10m) ## [1] 289.6244# df BIC
# PAR_LAI_mean_10m 7 290.8233
# mvPAR_LAI_mean_10m 7 289.6244
####PAR LAI sd ####
#p-value: 0.02881
PAR_LAI_sd_10m<-lm(PAR_LAI_sd~ eastwest,
data=bighornData)
summary(PAR_LAI_sd_10m)##
## Call:
## lm(formula = PAR_LAI_sd ~ eastwest, data = bighornData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73833 -0.27252 -0.06329 0.20379 1.57294
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.64220 0.04844 13.259 <2e-16 ***
## eastwest -0.10883 0.04877 -2.231 0.0288 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4138 on 71 degrees of freedom
## Multiple R-squared: 0.06553, Adjusted R-squared: 0.05237
## F-statistic: 4.979 on 1 and 71 DF, p-value: 0.02881AIC(PAR_LAI_sd_10m)## [1] 82.32141BIC(PAR_LAI_sd_10m) ## [1] 89.19279# mvPAR_LAI_sd_10m<-glm(PAR_LAI_sd~ eastwest,
# family=gaussian(link="log"),
# data=bighornData)
# summary(mvPAR_LAI_sd10m)
# BIC(PAR_LAI_sd_10m,mvPAR_LAI_sd_10m)
# df BIC
# PAR_LAI_sd_10m 3 89.19279
# mvPAR_LAI_sd_10m 3 89.71298
#####################################################################################
####PAR Average mean ####
#p-value: 0.008254
PAR_Average_mean_10m<-lm(PAR_Average_mean~ lulc_reduced,
data=bighornData)
summary(PAR_Average_mean_10m)##
## Call:
## lm(formula = PAR_Average_mean ~ lulc_reduced, data = bighornData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -410.28 -127.40 -67.51 107.00 777.92
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 420.48 61.84 6.799 3.06e-09 ***
## lulc_reduced149 -227.87 69.76 -3.266 0.0017 **
## lulc_reduced151 -286.57 94.46 -3.034 0.0034 **
## lulc_reduced155 -219.35 97.78 -2.243 0.0281 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 214.2 on 69 degrees of freedom
## Multiple R-squared: 0.1557, Adjusted R-squared: 0.119
## F-statistic: 4.241 on 3 and 69 DF, p-value: 0.008254densityPlot(PAR_Average_mean_10m$residuals)
AIC(PAR_Average_mean_10m)## [1] 996.6328BIC(PAR_Average_mean_10m) ## [1] 1008.085# mvPAR_Average_mean_10m<-glm(PAR_Average_mean~ lulc_reduced,
# family=gaussian(link="log"),
# data=bighornData)
# summary(mvPAR_Average_mean_10m)
# BIC(PAR_Average_mean_10m,mvPAR_Average_mean_10m)
# df BIC
# PAR_Average_mean_10m 5 1008.085
# mvPAR_Average_mean_10m 5 1008.085
####PAR Average SD ####
#p-value: 0.01623
# PAR_Average_sd10m<-lm(PAR_Average_sd~ TRI+slope+ northsouth,
# data=bighornData)
# summary(PAR_Average_sd10m)
mvPAR_Average_sd10m<-glm(PAR_Average_sd~ TRI+slope+ northsouth,
family=gaussian(link="log"),
data=bighornData)
summary(mvPAR_Average_sd10m)##
## Call:
## glm(formula = PAR_Average_sd ~ TRI + slope + northsouth, family = gaussian(link = "log"),
## data = bighornData)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -262.81 -104.18 -8.98 96.41 539.49
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.7304 0.2414 19.593 <2e-16 ***
## TRI 1.6649 0.8282 2.010 0.0483 *
## slope -2.2299 0.9365 -2.381 0.0200 *
## northsouth -0.3935 0.1950 -2.018 0.0475 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 28999.74)
##
## Null deviance: 2451433 on 72 degrees of freedom
## Residual deviance: 2000577 on 69 degrees of freedom
## AIC: 963.11
##
## Number of Fisher Scoring iterations: 17AIC(mvPAR_Average_sd10m)## [1] 963.1146BIC(mvPAR_Average_sd10m) ## [1] 974.5668# BIC(PAR_Average_sd10m,mvPAR_Average_sd10m)
# df BIC
# PAR_Average_sd10m 5 978.5905
# mvPAR_Average_sd10m 5 974.5668
#####################################################################################
#### Canopy Density Mean ####
#p-value: 8.287e-06
Canopy_Density_mean_10m<-lm(Canopy_Density_mean~ DEM+ lulc_reduced,
data=bighornData)
summary(Canopy_Density_mean_10m)##
## Call:
## lm(formula = Canopy_Density_mean ~ DEM + lulc_reduced, data = bighornData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.178 -2.206 1.187 4.191 22.762
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 55.305 2.925 18.908 < 2e-16 ***
## DEM 3.612 1.328 2.721 0.00826 **
## lulc_reduced149 14.206 3.317 4.283 5.93e-05 ***
## lulc_reduced151 12.912 4.605 2.804 0.00658 **
## lulc_reduced155 8.106 4.715 1.719 0.09010 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.06 on 68 degrees of freedom
## Multiple R-squared: 0.3421, Adjusted R-squared: 0.3034
## F-statistic: 8.841 on 4 and 68 DF, p-value: 8.287e-06AIC(Canopy_Density_mean_10m)## [1] 551.0652BIC(Canopy_Density_mean_10m)## [1] 564.8079densityPlot(Canopy_Density_mean_10m$residuals)
# mvCanopy_Density_mean_10m<-glm(Canopy_Density_mean~ DEM+ lulc_reduced,
# family=gaussian(link="logit"),
# data=bighornData)
# summary(mvCanopy_Density_mean_10m)
#BIC(Canopy_Density_mean_10m,mvCanopy_Density_mean_10m)
# df BIC
# Canopy_Density_mean_10m 6 564.8079
# mvCanopy_Density_mean_10m 6 566.1761
#####Canopy Density SD ####
#5.175e-05
# Canopy_Density_sd_10m<-lm(Canopy_Density_sd~ DEM+tpi100m+ lulc_reduced,
# data=bighornData)
# summary(Canopy_Density_sd_10m)
mvCanopy_Density_sd_10m<-glm(Canopy_Density_sd~ DEM+tpi100m+ lulc_reduced,
family=gaussian(link="log"),
data=bighornData)
summary(mvCanopy_Density_sd_10m)##
## Call:
## glm(formula = Canopy_Density_sd ~ DEM + tpi100m + lulc_reduced,
## family = gaussian(link = "log"), data = bighornData)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.4643 -2.3473 -1.1044 0.4404 14.6721
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.23457 0.15057 14.841 < 2e-16 ***
## DEM -0.16665 0.07603 -2.192 0.03186 *
## tpi100m -0.25758 0.09212 -2.796 0.00674 **
## lulc_reduced149 -0.74751 0.22374 -3.341 0.00137 **
## lulc_reduced151 -0.88263 0.45086 -1.958 0.05444 .
## lulc_reduced155 -0.61829 0.28028 -2.206 0.03082 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 22.61968)
##
## Null deviance: 2338.1 on 72 degrees of freedom
## Residual deviance: 1515.5 on 67 degrees of freedom
## AIC: 442.58
##
## Number of Fisher Scoring iterations: 8AIC(mvCanopy_Density_sd_10m)## [1] 442.5765BIC(mvCanopy_Density_sd_10m)## [1] 458.6098#BIC(Canopy_Density_sd_10m,mvCanopy_Density_sd_10m)
# df BIC
# Canopy_Density_sd_10m 7 461.2300
# mvCanopy_Density_sd_10m 7 458.6098
####################################################################################
#### Leaf Angle mean ####
#p-value: 0.5935
Leaf_Angle_mean_10<-lm(Leaf_Angle_mean~ lulc_reduced,
data=bighornData)
summary(Leaf_Angle_mean_10)##
## Call:
## lm(formula = Leaf_Angle_mean ~ lulc_reduced, data = bighornData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.627 -8.346 1.694 9.114 35.162
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.727 4.390 11.100 <2e-16 ***
## lulc_reduced149 5.660 4.952 1.143 0.257
## lulc_reduced151 1.049 6.706 0.156 0.876
## lulc_reduced155 1.311 6.941 0.189 0.851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.21 on 69 degrees of freedom
## Multiple R-squared: 0.02696, Adjusted R-squared: -0.01534
## F-statistic: 0.6374 on 3 and 69 DF, p-value: 0.5935AIC(Leaf_Angle_mean_10)## [1] 610.4254BIC(Leaf_Angle_mean_10)## [1] 621.8777densityPlot(Leaf_Angle_mean_10$residuals)
#### Leaf Angle sd ####
#p-value: 0.001394
Leaf_Angle_sd_10<-lm(Leaf_Angle_sd~ tpi1000m+lulc_reduced,
data=bighornData)
summary(Leaf_Angle_sd_10)##
## Call:
## lm(formula = Leaf_Angle_sd ~ tpi1000m + lulc_reduced, data = bighornData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.8374 -6.5164 -0.5118 7.8988 24.1936
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.7175 2.9241 8.453 3.67e-12 ***
## tpi1000m -4.1758 1.2449 -3.354 0.00131 **
## lulc_reduced149 -6.5279 3.3152 -1.969 0.05308 .
## lulc_reduced151 -0.9586 4.4813 -0.214 0.83126
## lulc_reduced155 -3.1830 4.6759 -0.681 0.49839
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.13 on 67 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2296, Adjusted R-squared: 0.1836
## F-statistic: 4.991 on 4 and 67 DF, p-value: 0.001394AIC(Leaf_Angle_sd_10)## [1] 544.5641BIC(Leaf_Angle_sd_10)## [1] 558.2241densityPlot(Leaf_Angle_sd_10$residuals)
#####################################################################################
####Diversity ####
# p-value: 0.0009821
Diversity_10m<-lm(div~DEM+northsouth+lulc_reduced,
data=bighornData)
summary(Diversity_10m) ##
## Call:
## lm(formula = div ~ DEM + northsouth + lulc_reduced, data = bighornData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3811 -0.1322 -0.0067 0.1049 0.7027
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.38439 0.06148 6.252 3.22e-08 ***
## DEM 0.09229 0.02705 3.412 0.0011 **
## northsouth -0.04444 0.02512 -1.769 0.0815 .
## lulc_reduced149 0.16987 0.06971 2.437 0.0175 *
## lulc_reduced151 0.03923 0.09604 0.408 0.6842
## lulc_reduced155 0.15174 0.09824 1.545 0.1272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.205 on 67 degrees of freedom
## Multiple R-squared: 0.2594, Adjusted R-squared: 0.2041
## F-statistic: 4.694 on 5 and 67 DF, p-value: 0.0009821AIC(Diversity_10m)## [1] -16.45697BIC(Diversity_10m)## [1] -0.4237531densityPlot(Diversity_10m$residuals)
#####################################################################################
#### Average DBH ####
#p-value: 0.0009843
# DBH_mean_10m<-lm(Average_DBH~ tpi100m+ lulc_reduced,
# data=bighornData)
# summary(DBH_mean_10m)
mvDBH_mean_10m<-glm(Average_DBH~ tpi100m+ lulc_reduced,
family=gaussian(link="log"),
data=bighornData)
summary(mvDBH_mean_10m)##
## Call:
## glm(formula = Average_DBH ~ tpi100m + lulc_reduced, family = gaussian(link = "log"),
## data = bighornData)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -18.515 -5.417 -1.854 5.174 28.218
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.44753 0.08187 42.110 < 2e-16 ***
## tpi100m -0.11882 0.03845 -3.091 0.00289 **
## lulc_reduced149 -0.25190 0.09853 -2.557 0.01281 *
## lulc_reduced151 0.10019 0.11714 0.855 0.39542
## lulc_reduced155 -0.29661 0.15219 -1.949 0.05542 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 77.95788)
##
## Null deviance: 7049.7 on 72 degrees of freedom
## Residual deviance: 5301.1 on 68 degrees of freedom
## AIC: 531.99
##
## Number of Fisher Scoring iterations: 5AIC(mvDBH_mean_10m)## [1] 531.9853BIC(mvDBH_mean_10m)## [1] 545.7281#BIC(DBH_mean_10m,mvDBH_mean_10m)
# df BIC
# DBH_mean_10m 6 546.9549
# mvDBH_mean_10m 6 545.7281
#####################################################################################
#### SD DBH ####
#p-value: 0.001216
# DBH_sd_10m<-lm(DBH_Sd~ tpi1000m+ lulc_reduced,
# data=bighornData)
# summary(DBH_sd_10m)
mvDBH_sd_10m<-glm(DBH_Sd~ tpi1000m+ lulc_reduced,
family=gaussian(link="log"),
start = c(1,1,1,1,1),
data=bighornData)
summary(mvDBH_sd_10m)##
## Call:
## glm(formula = DBH_Sd ~ tpi1000m + lulc_reduced, family = gaussian(link = "log"),
## data = bighornData, start = c(1, 1, 1, 1, 1))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -10.9495 -2.9485 -0.4068 3.0472 15.9326
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.91597 0.22261 8.607 1.94e-12 ***
## tpi1000m -0.20102 0.07326 -2.744 0.00778 **
## lulc_reduced149 0.14843 0.24541 0.605 0.54732
## lulc_reduced151 0.76101 0.24972 3.047 0.00330 **
## lulc_reduced155 -0.22064 0.37252 -0.592 0.55565
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 29.08617)
##
## Null deviance: 2658.4 on 71 degrees of freedom
## Residual deviance: 1948.7 on 67 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 453.8
##
## Number of Fisher Scoring iterations: 10AIC(mvDBH_sd_10m)## [1] 453.803BIC(mvDBH_sd_10m)## [1] 467.463densityPlot(mvDBH_sd_10m$residuals)
#BIC(DBH_sd_10m,mvDBH_sd_10m)
# df BIC
# DBH_sd_10m 6 470.7204
# mvDBH_sd_10m 6 467.4630
#####################################################################################
#S####tand Density ####
#p-value: 1.806e-06
Stand_Density_10m<-lm(Stand_Density~ DEM+ tpi100m+ lulc_reduced,
data=bighornData)
summary(Stand_Density_10m) ##
## Call:
## lm(formula = Stand_Density ~ DEM + tpi100m + lulc_reduced, data = bighornData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9979 -1.8267 -0.1134 1.6456 9.6533
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.0015 0.9799 6.125 5.38e-08 ***
## DEM 1.2977 0.4563 2.844 0.00591 **
## tpi100m 1.3722 0.4234 3.241 0.00186 **
## lulc_reduced149 3.0350 1.1082 2.739 0.00790 **
## lulc_reduced151 -0.2019 1.5491 -0.130 0.89671
## lulc_reduced155 1.5211 1.5988 0.951 0.34483
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.36 on 67 degrees of freedom
## Multiple R-squared: 0.3969, Adjusted R-squared: 0.3519
## F-statistic: 8.819 on 5 and 67 DF, p-value: 1.806e-06AIC(Stand_Density_10m)## [1] 391.8308BIC(Stand_Density_10m)## [1] 407.864densityPlot(Stand_Density_10m$residuals)
#####################################################################################Tree Height Regressions at 10 Meters
#Mean Tree Height
#p-value: 0.0199
Mean_Tree_Height10m<-lm(Mean_Tree_Height~northsouth+lulc_reduced,
data=bighornDataTreeHeightsAndAge)
summary(Mean_Tree_Height10m)##
## Call:
## lm(formula = Mean_Tree_Height ~ northsouth + lulc_reduced, data = bighornDataTreeHeightsAndAge)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.9490 -2.0172 -0.4666 2.3026 9.2168
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.3446 1.3570 5.412 3.64e-06 ***
## northsouth -0.8025 0.5636 -1.424 0.16268
## lulc_reduced149 5.3506 1.5301 3.497 0.00122 **
## lulc_reduced151 5.5459 2.3165 2.394 0.02170 *
## lulc_reduced155 6.9510 3.6681 1.895 0.06572 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.43 on 38 degrees of freedom
## (30 observations deleted due to missingness)
## Multiple R-squared: 0.259, Adjusted R-squared: 0.181
## F-statistic: 3.321 on 4 and 38 DF, p-value: 0.0199AIC(Mean_Tree_Height10m)## [1] 234.7079BIC(Mean_Tree_Height10m)## [1] 245.2751densityPlot( Mean_Tree_Height10m$residuals)
###################################################################################
#SD Tree Height
#p-value: 0.006275
SD_Tree_Height10m<-lm(Sd_Tree_Height~ tpi100m+lulc_reduced,
data=bighornDataTreeHeightsAndAge)
summary(SD_Tree_Height10m)##
## Call:
## lm(formula = Sd_Tree_Height ~ tpi100m + lulc_reduced, data = bighornDataTreeHeightsAndAge)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0817 -1.2138 -0.1103 0.5723 4.4815
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0266 0.6550 7.674 3.07e-09 ***
## tpi100m -0.8031 0.2634 -3.049 0.00417 **
## lulc_reduced149 -1.0080 0.7271 -1.386 0.17371
## lulc_reduced151 0.7790 1.0901 0.715 0.47922
## lulc_reduced155 -0.8551 1.8526 -0.462 0.64704
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.732 on 38 degrees of freedom
## (30 observations deleted due to missingness)
## Multiple R-squared: 0.308, Adjusted R-squared: 0.2352
## F-statistic: 4.229 on 4 and 38 DF, p-value: 0.006275AIC(SD_Tree_Height10m)## [1] 175.9737BIC(SD_Tree_Height10m)## [1] 186.5409densityPlot(SD_Tree_Height10m$residuals)
##################################################################################
#Max Tree Height (northsouth aspect)
#p-value: 0.00969
Max_Tree_Height10m<-lm(Max_Tree_Height~ aspect+ tpi1000m +lulc_reduced,
data=bighornDataTreeHeightsAndAge)
summary(Max_Tree_Height10m)##
## Call:
## lm(formula = Max_Tree_Height ~ aspect + tpi1000m + lulc_reduced,
## data = bighornDataTreeHeightsAndAge)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.6333 -2.8068 0.3041 1.9522 6.8380
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.8762 1.3137 12.085 2.07e-14 ***
## aspect 0.8112 0.5026 1.614 0.1150
## tpi1000m -1.1283 0.5897 -1.913 0.0635 .
## lulc_reduced149 2.3294 1.4470 1.610 0.1159
## lulc_reduced151 5.5133 2.1752 2.535 0.0156 *
## lulc_reduced155 6.5971 3.6305 1.817 0.0773 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.388 on 37 degrees of freedom
## (30 observations deleted due to missingness)
## Multiple R-squared: 0.326, Adjusted R-squared: 0.235
## F-statistic: 3.58 on 5 and 37 DF, p-value: 0.00969AIC(Max_Tree_Height10m)## [1] 234.5193BIC(Max_Tree_Height10m)## [1] 246.8477densityPlot(Max_Tree_Height10m$residuals)
##################################################################################
#Average Stand Age
#p-value: 0.0001305
Average_Stand_Age10m<-lm(Average_Stand_Age~aspect+ tpi1000m +lulc_reduced,
data=bighornDataTreeHeightsAndAge)
summary(Average_Stand_Age10m)##
## Call:
## lm(formula = Average_Stand_Age ~ aspect + tpi1000m + lulc_reduced,
## data = bighornDataTreeHeightsAndAge)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.865 -10.901 -2.243 8.560 56.634
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 74.491 5.781 12.885 < 2e-16 ***
## aspect 5.211 2.422 2.151 0.03510 *
## tpi1000m -6.481 2.449 -2.646 0.01016 *
## lulc_reduced149 -19.671 6.534 -3.010 0.00369 **
## lulc_reduced151 1.676 8.864 0.189 0.85066
## lulc_reduced155 -18.959 9.432 -2.010 0.04852 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.89 on 66 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.3113, Adjusted R-squared: 0.2591
## F-statistic: 5.965 on 5 and 66 DF, p-value: 0.0001305AIC(Average_Stand_Age10m)## [1] 642.6376BIC(Average_Stand_Age10m)## [1] 658.5743densityPlot(Average_Stand_Age10m$residuals)
Regression Results
Regression Results
#Gap Fraction Mean
par(mfrow=c(2,3))
plot(GapFraction_mean_10m, 1:6)
gvlma(GapFraction_mean_10m)##
## Call:
## lm(formula = Gap_Fraction_mean ~ TRI + tpi1000m + lulc_reduced,
## data = bighornData)
##
## Coefficients:
## (Intercept) TRI tpi1000m lulc_reduced149
## 1.11244 0.07958 0.08333 0.49340
## lulc_reduced151 lulc_reduced155
## 0.60939 0.16545
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = GapFraction_mean_10m)
##
## Value p-value Decision
## Global Stat 12.885 0.01185 Assumptions NOT satisfied!
## Skewness 1.457 0.22742 Assumptions acceptable.
## Kurtosis 1.715 0.19028 Assumptions acceptable.
## Link Function 4.881 0.02716 Assumptions NOT satisfied!
## Heteroscedasticity 4.832 0.02794 Assumptions NOT satisfied!# Coefficients:
# (Intercept) TRI tpi1000m lulc_reduced149 lulc_reduced151 lulc_reduced155
# 1.11244 0.07958 0.08333 0.49340 0.60939 0.16545
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = GapFraction_mean_10m)
#
# Value p-value Decision
# Global Stat 12.885 0.01185 Assumptions NOT satisfied!
# Skewness 1.457 0.22742 Assumptions acceptable.
# Kurtosis 1.715 0.19028 Assumptions acceptable.
# Link Function 4.881 0.02716 Assumptions NOT satisfied!
# Heteroscedasticity 4.832 0.02794 Assumptions NOT satisfied!
#Gap Fraction SD
par(mfrow=c(2,3))
plot(MVGapFraction_sd_10m, 1:6)
#PAR LAI Mean
par(mfrow=c(2,3))
plot(mvPAR_LAI_mean_10m, 1:6)
#PAR LAI SD
par(mfrow=c(2,3))
plot(PAR_LAI_sd_10m, 1:6) 
gvlma(PAR_LAI_sd_10m)##
## Call:
## lm(formula = PAR_LAI_sd ~ eastwest, data = bighornData)
##
## Coefficients:
## (Intercept) eastwest
## 0.6422 -0.1088
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = PAR_LAI_sd_10m)
##
## Value p-value Decision
## Global Stat 31.2255 2.754e-06 Assumptions NOT satisfied!
## Skewness 15.6306 7.700e-05 Assumptions NOT satisfied!
## Kurtosis 14.2828 1.573e-04 Assumptions NOT satisfied!
## Link Function 0.7860 3.753e-01 Assumptions acceptable.
## Heteroscedasticity 0.5261 4.683e-01 Assumptions acceptable.# Call:
# lm(formula = PAR_LAI_sd ~ eastwest, data = bighornData)
#
# Coefficients:
# (Intercept) eastwest
# 0.6422 -0.1088
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = PAR_LAI_sd_10m)
#
# Value p-value Decision
# Global Stat 31.2255 2.754e-06 Assumptions NOT satisfied!
# Skewness 15.6306 7.700e-05 Assumptions NOT satisfied!
# Kurtosis 14.2828 1.573e-04 Assumptions NOT satisfied!
# Link Function 0.7860 3.753e-01 Assumptions acceptable.
# Heteroscedasticity 0.5261 4.683e-01 Assumptions acceptable.
#PAR Average Mean
par(mfrow=c(2,3))
plot(PAR_Average_mean_10m, 1:6)
gvlma(PAR_Average_mean_10m)##
## Call:
## lm(formula = PAR_Average_mean ~ lulc_reduced, data = bighornData)
##
## Coefficients:
## (Intercept) lulc_reduced149 lulc_reduced151 lulc_reduced155
## 420.5 -227.9 -286.6 -219.4
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = PAR_Average_mean_10m)
##
## Value p-value Decision
## Global Stat 4.265e+01 1.222e-08 Assumptions NOT satisfied!
## Skewness 1.811e+01 2.089e-05 Assumptions NOT satisfied!
## Kurtosis 1.321e+01 2.779e-04 Assumptions NOT satisfied!
## Link Function 1.853e-14 1.000e+00 Assumptions acceptable.
## Heteroscedasticity 1.133e+01 7.623e-04 Assumptions NOT satisfied!# Call:
# lm(formula = PAR_Average_mean ~ lulc_reduced, data = bighornData)
#
# Coefficients:
# (Intercept) lulc_reduced149 lulc_reduced151 lulc_reduced155
# 420.5 -227.9 -286.6 -219.4
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = PAR_Average_mean_10m)
#
# Value p-value Decision
# Global Stat 4.265e+01 1.222e-08 Assumptions NOT satisfied!
# Skewness 1.811e+01 2.089e-05 Assumptions NOT satisfied!
# Kurtosis 1.321e+01 2.779e-04 Assumptions NOT satisfied!
# Link Function 1.853e-14 1.000e+00 Assumptions acceptable.
# Heteroscedasticity 1.133e+01 7.623e-04 Assumptions NOT satisfied!
#PAR Average SD
par(mfrow=c(2,3))
plot(mvPAR_Average_sd10m, 1:6)
#Canopy Density Mean
par(mfrow=c(2,3))
plot(Canopy_Density_mean_10m, 1:6)
gvlma(Canopy_Density_mean_10m)##
## Call:
## lm(formula = Canopy_Density_mean ~ DEM + lulc_reduced, data = bighornData)
##
## Coefficients:
## (Intercept) DEM lulc_reduced149 lulc_reduced151
## 55.305 3.612 14.206 12.912
## lulc_reduced155
## 8.106
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = Canopy_Density_mean_10m)
##
## Value p-value Decision
## Global Stat 87.465 0.000e+00 Assumptions NOT satisfied!
## Skewness 23.079 1.555e-06 Assumptions NOT satisfied!
## Kurtosis 51.359 7.694e-13 Assumptions NOT satisfied!
## Link Function 4.812 2.826e-02 Assumptions NOT satisfied!
## Heteroscedasticity 8.216 4.153e-03 Assumptions NOT satisfied!# Call:
# lm(formula = Canopy_Density_mean ~ DEM + lulc_reduced, data = bighornData)
#
# Coefficients:
# (Intercept) DEM lulc_reduced149 lulc_reduced151 lulc_reduced155
# 55.305 3.612 14.206 12.912 8.106
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = Canopy_Density_mean_10m)
#
# Value p-value Decision
# Global Stat 87.465 0.000e+00 Assumptions NOT satisfied!
# Skewness 23.079 1.555e-06 Assumptions NOT satisfied!
# Kurtosis 51.359 7.694e-13 Assumptions NOT satisfied!
# Link Function 4.812 2.826e-02 Assumptions NOT satisfied!
# Heteroscedasticity 8.216 4.153e-03 Assumptions NOT satisfied!
#Canopy Density SD
par(mfrow=c(2,3))
plot(mvCanopy_Density_sd_10m, 1:6)
#Leaf Angle Mean
par(mfrow=c(2,3))
plot(Leaf_Angle_mean_10, 1:6)
gvlma(Leaf_Angle_mean_10)##
## Call:
## lm(formula = Leaf_Angle_mean ~ lulc_reduced, data = bighornData)
##
## Coefficients:
## (Intercept) lulc_reduced149 lulc_reduced151 lulc_reduced155
## 48.727 5.660 1.049 1.311
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = Leaf_Angle_mean_10)
##
## Value p-value Decision
## Global Stat 6.230e+00 0.1826 Assumptions acceptable.
## Skewness 2.686e+00 0.1012 Assumptions acceptable.
## Kurtosis 1.936e+00 0.1641 Assumptions acceptable.
## Link Function 5.268e-13 1.0000 Assumptions acceptable.
## Heteroscedasticity 1.608e+00 0.2048 Assumptions acceptable.# Call:
# lm(formula = Leaf_Angle_mean ~ lulc_reduced, data = bighornData)
#
# Coefficients:
# (Intercept) lulc_reduced149 lulc_reduced151 lulc_reduced155
# 48.727 5.660 1.049 1.311
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = Leaf_Angle_mean_10)
#
# Value p-value Decision
# Global Stat 6.230e+00 0.1826 Assumptions acceptable.
# Skewness 2.686e+00 0.1012 Assumptions acceptable.
# Kurtosis 1.936e+00 0.1641 Assumptions acceptable.
# Link Function 5.268e-13 1.0000 Assumptions acceptable.
# Heteroscedasticity 1.608e+00 0.2048 Assumptions acceptable.
#Leaf Angle SD
par(mfrow=c(2,3))
plot(Leaf_Angle_sd_10, 1:6)
gvlma(Leaf_Angle_sd_10)##
## Call:
## lm(formula = Leaf_Angle_sd ~ tpi1000m + lulc_reduced, data = bighornData)
##
## Coefficients:
## (Intercept) tpi1000m lulc_reduced149 lulc_reduced151
## 24.7175 -4.1758 -6.5279 -0.9586
## lulc_reduced155
## -3.1830
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = Leaf_Angle_sd_10)
##
## Value p-value Decision
## Global Stat 4.89145 0.2986 Assumptions acceptable.
## Skewness 0.48586 0.4858 Assumptions acceptable.
## Kurtosis 0.71337 0.3983 Assumptions acceptable.
## Link Function 3.67004 0.0554 Assumptions acceptable.
## Heteroscedasticity 0.02219 0.8816 Assumptions acceptable.# Call:
# lm(formula = Leaf_Angle_sd ~ tpi1000m + lulc_reduced, data = bighornData)
#
# Coefficients:
# (Intercept) tpi1000m lulc_reduced149 lulc_reduced151 lulc_reduced155
# 24.7175 -4.1758 -6.5279 -0.9586 -3.1830
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = Leaf_Angle_sd_10)
#
# Value p-value Decision
# Global Stat 4.89145 0.2986 Assumptions acceptable.
# Skewness 0.48586 0.4858 Assumptions acceptable.
# Kurtosis 0.71337 0.3983 Assumptions acceptable.
# Link Function 3.67004 0.0554 Assumptions acceptable.
# Heteroscedasticity 0.02219 0.8816 Assumptions acceptable
#Diversity
par(mfrow=c(2,3))
plot(Diversity_10m, 1:6)
gvlma(Diversity_10m)##
## Call:
## lm(formula = div ~ DEM + northsouth + lulc_reduced, data = bighornData)
##
## Coefficients:
## (Intercept) DEM northsouth lulc_reduced149
## 0.38439 0.09229 -0.04444 0.16987
## lulc_reduced151 lulc_reduced155
## 0.03923 0.15174
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = Diversity_10m)
##
## Value p-value Decision
## Global Stat 17.0809 0.001864 Assumptions NOT satisfied!
## Skewness 6.2070 0.012725 Assumptions NOT satisfied!
## Kurtosis 5.4838 0.019194 Assumptions NOT satisfied!
## Link Function 0.8535 0.355558 Assumptions acceptable.
## Heteroscedasticity 4.5366 0.033178 Assumptions NOT satisfied!# Call:
# lm(formula = div ~ DEM + northsouth + lulc_reduced, data = bighornData)
#
# Coefficients:
# (Intercept) DEM northsouth lulc_reduced149 lulc_reduced151
# 0.38439 0.09229 -0.04444 0.16987 0.03923
# lulc_reduced155
# 0.15174
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = Diversity_10m)
#
# Value p-value Decision
# Global Stat 17.0809 0.001864 Assumptions NOT satisfied!
# Skewness 6.2070 0.012725 Assumptions NOT satisfied!
# Kurtosis 5.4838 0.019194 Assumptions NOT satisfied!
# Link Function 0.8535 0.355558 Assumptions acceptable.
# Heteroscedasticity 4.5366 0.033178 Assumptions NOT satisfied!
#DBH Mean
par(mfrow=c(2,3))
plot(mvDBH_mean_10m, 1:6)
#DBH SD
par(mfrow=c(2,3))
plot(mvDBH_sd_10m, 1:6)
#Stand Density
par(mfrow=c(2,3))
plot(Stand_Density_10m, 1:6)
gvlma(Stand_Density_10m)##
## Call:
## lm(formula = Stand_Density ~ DEM + tpi100m + lulc_reduced, data = bighornData)
##
## Coefficients:
## (Intercept) DEM tpi100m lulc_reduced149
## 6.0015 1.2977 1.3722 3.0350
## lulc_reduced151 lulc_reduced155
## -0.2019 1.5211
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = Stand_Density_10m)
##
## Value p-value Decision
## Global Stat 8.1492 0.08626 Assumptions acceptable.
## Skewness 2.6182 0.10564 Assumptions acceptable.
## Kurtosis 1.3827 0.23964 Assumptions acceptable.
## Link Function 3.3454 0.06739 Assumptions acceptable.
## Heteroscedasticity 0.8029 0.37023 Assumptions acceptable.# Call:
# lm(formula = Stand_Density ~ DEM + tpi100m + lulc_reduced, data = bighornData)
#
# Coefficients:
# (Intercept) DEM tpi100m lulc_reduced149 lulc_reduced151 lulc_reduced155
# 6.0015 1.2977 1.3722 3.0350 -0.2019 1.5211
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = Stand_Density_10m)
#
# Value p-value Decision
# Global Stat 8.1492 0.08626 Assumptions acceptable.
# Skewness 2.6182 0.10564 Assumptions acceptable.
# Kurtosis 1.3827 0.23964 Assumptions acceptable.
# Link Function 3.3454 0.06739 Assumptions acceptable.
# Heteroscedasticity 0.8029 0.37023 Assumptions acceptable.
# Tree Height Mean
par(mfrow=c(2,3))
plot(Mean_Tree_Height10m, 1:6)
gvlma(Mean_Tree_Height10m)##
## Call:
## lm(formula = Mean_Tree_Height ~ northsouth + lulc_reduced, data = bighornDataTreeHeightsAndAge)
##
## Coefficients:
## (Intercept) northsouth lulc_reduced149 lulc_reduced151
## 7.3446 -0.8025 5.3506 5.5459
## lulc_reduced155
## 6.9510
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = Mean_Tree_Height10m)
##
## Value p-value Decision
## Global Stat 1.052881 0.9017 Assumptions acceptable.
## Skewness 0.997379 0.3179 Assumptions acceptable.
## Kurtosis 0.036870 0.8477 Assumptions acceptable.
## Link Function 0.005801 0.9393 Assumptions acceptable.
## Heteroscedasticity 0.012832 0.9098 Assumptions acceptable.# Call:
# lm(formula = Mean_Tree_Height ~ northsouth + lulc_reduced, data = bighornDataTreeHeightsAndAge)
#
# Coefficients:
# (Intercept) northsouth lulc_reduced149 lulc_reduced151 lulc_reduced155
# 7.3446 -0.8025 5.3506 5.5459 6.9510
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = Mean_Tree_Height10m)
#
# Value p-value Decision
# Global Stat 1.052881 0.9017 Assumptions acceptable.
# Skewness 0.997379 0.3179 Assumptions acceptable.
# Kurtosis 0.036870 0.8477 Assumptions acceptable.
# Link Function 0.005801 0.9393 Assumptions acceptable.
# Heteroscedasticity 0.012832 0.9098 Assumptions acceptable.
# Tree Height SD
par(mfrow=c(2,3))
plot(SD_Tree_Height10m, 1:6)
gvlma(SD_Tree_Height10m)##
## Call:
## lm(formula = Sd_Tree_Height ~ tpi100m + lulc_reduced, data = bighornDataTreeHeightsAndAge)
##
## Coefficients:
## (Intercept) tpi100m lulc_reduced149 lulc_reduced151
## 5.0266 -0.8031 -1.0080 0.7790
## lulc_reduced155
## -0.8551
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = SD_Tree_Height10m)
##
## Value p-value Decision
## Global Stat 4.9753 0.28984 Assumptions acceptable.
## Skewness 3.9299 0.04744 Assumptions NOT satisfied!
## Kurtosis 0.3955 0.52942 Assumptions acceptable.
## Link Function 0.3909 0.53182 Assumptions acceptable.
## Heteroscedasticity 0.2590 0.61080 Assumptions acceptable.# Call:
# lm(formula = Sd_Tree_Height ~ tpi100m + lulc_reduced, data = bighornDataTreeHeightsAndAge)
#
# Coefficients:
# (Intercept) tpi100m lulc_reduced149 lulc_reduced151 lulc_reduced155
# 5.0266 -0.8031 -1.0080 0.7790 -0.8551
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = SD_Tree_Height10m)
#
# Value p-value Decision
# Global Stat 4.9753 0.28984 Assumptions acceptable.
# Skewness 3.9299 0.04744 Assumptions NOT satisfied!
# Kurtosis 0.3955 0.52942 Assumptions acceptable.
# Link Function 0.3909 0.53182 Assumptions acceptable.
# Heteroscedasticity 0.2590 0.61080 Assumptions acceptable.
#Max Tree Height
par(mfrow=c(2,3))
plot(Max_Tree_Height10m, 1:6)
gvlma(Max_Tree_Height10m)##
## Call:
## lm(formula = Max_Tree_Height ~ aspect + tpi1000m + lulc_reduced,
## data = bighornDataTreeHeightsAndAge)
##
## Coefficients:
## (Intercept) aspect tpi1000m lulc_reduced149
## 15.8762 0.8112 -1.1283 2.3294
## lulc_reduced151 lulc_reduced155
## 5.5133 6.5971
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = Max_Tree_Height10m)
##
## Value p-value Decision
## Global Stat 4.26114 0.37182 Assumptions acceptable.
## Skewness 0.46840 0.49372 Assumptions acceptable.
## Kurtosis 0.85119 0.35622 Assumptions acceptable.
## Link Function 2.91962 0.08751 Assumptions acceptable.
## Heteroscedasticity 0.02193 0.88228 Assumptions acceptable.# Call:
# lm(formula = Max_Tree_Height ~ aspect + tpi1000m + lulc_reduced,
# data = bighornDataTreeHeightsAndAge)
#
# Coefficients:
# (Intercept) aspect tpi1000m lulc_reduced149 lulc_reduced151 lulc_reduced155
# 15.8762 0.8112 -1.1283 2.3294 5.5133 6.5971
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = Max_Tree_Height10m)
#
# Value p-value Decision
# Global Stat 4.26114 0.37182 Assumptions acceptable.
# Skewness 0.46840 0.49372 Assumptions acceptable.
# Kurtosis 0.85119 0.35622 Assumptions acceptable.
# Link Function 2.91962 0.08751 Assumptions acceptable.
# Heteroscedasticity 0.02193 0.88228 Assumptions acceptable.
#Average Stand Age
par(mfrow=c(2,3))
plot(Average_Stand_Age10m, 1:6)
gvlma(Average_Stand_Age10m)##
## Call:
## lm(formula = Average_Stand_Age ~ aspect + tpi1000m + lulc_reduced,
## data = bighornDataTreeHeightsAndAge)
##
## Coefficients:
## (Intercept) aspect tpi1000m lulc_reduced149
## 74.491 5.211 -6.481 -19.671
## lulc_reduced151 lulc_reduced155
## 1.676 -18.959
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = Average_Stand_Age10m)
##
## Value p-value Decision
## Global Stat 15.3846 0.003966 Assumptions NOT satisfied!
## Skewness 4.6874 0.030385 Assumptions NOT satisfied!
## Kurtosis 1.2152 0.270299 Assumptions acceptable.
## Link Function 0.8314 0.361869 Assumptions acceptable.
## Heteroscedasticity 8.6506 0.003270 Assumptions NOT satisfied!# Call:
# lm(formula = Average_Stand_Age ~ aspect + tpi1000m + lulc_reduced,
# data = bighornDataTreeHeightsAndAge)
#
# Coefficients:
# (Intercept) aspect tpi1000m lulc_reduced149 lulc_reduced151 lulc_reduced155
# 74.491 5.211 -6.481 -19.671 l1.676 -18.959
#
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance = 0.05
#
# Call:
# gvlma(x = Average_Stand_Age10m)
#
# Value p-value Decision
# Global Stat 15.3846 0.003966 Assumptions NOT satisfied!
# Skewness 4.6874 0.030385 Assumptions NOT satisfied!
# Kurtosis 1.2152 0.270299 Assumptions acceptable.
# Link Function 0.8314 0.361869 Assumptions acceptable.
# Heteroscedasticity 8.6506 0.003270 Assumptions NOT satisfied!Gap Fraction Plots
#GapFraction_mean_10m<-lm(Gap_Fraction_mean~ TRI+ tpi1000m+lulc_reduced,data=bighornData)
GF.TRI<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=TRI, y=Gap_Fraction_mean, color=Gap_Fraction_mean)) +
scale_color_viridis_c(option="turbo")+
labs(x = "Topographic Roughness Index", y ="Gap Fraction Mean",
title ="Gap Fraction Mean vs Topographic Roughness IndexI", color ="Gap Fraction Mean") +
theme(plot.title = element_text(hjust = 0.5))
GF.TRI
GF.100TPI<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=tpi100m, y=Gap_Fraction_mean, color=Gap_Fraction_mean)) +
scale_color_viridis_c(option="turbo")+
labs(x = "100m TPI", y ="Gap Fraction Mean",
title ="Gap Fraction Mean vs 100m TPI", color ="Gap Fraction Mean") +
theme(plot.title = element_text(hjust = 0.5))
GF.100TPI
GF.LULC<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=Gap_Fraction_mean, color=Gap_Fraction_mean)) +
scale_color_viridis_c(option="turbo")+
labs(x = "Landcover Type", y ="Gap Fraction Mean",
title ="Gap Fraction Mean vs Landcover Type", color ="Gap Fraction Mean") +
theme(plot.title = element_text(hjust = 0.5))
GF.LULC
Gap_Fraction_mean_Car<-car::avPlots(model = GapFraction_mean_10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=TRUE,
main=("Gap Fraction Added-Variable Plots")) 
#GapFraction_sd_10m<-lm(Gap_Fraction_sd~ DEM, data=bighornData)
GFsd.DEM<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=DEM, y=Gap_Fraction_sd, color=Gap_Fraction_sd)) +
scale_color_viridis_c(option="turbo")+
labs(x = "Elevation", y ="Gap Fraction Sd",
title ="Gap Fraction Sd vs Elevation", color ="Gap Fraction Sd") +
theme(plot.title = element_text(hjust = 0.5))
GFsd.DEM
Gap_Fraction_sd_Car<-car::avPlots(model = MVGapFraction_sd_10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Gap Fraction Added-Variable Plot")) 
PAR LAI Plots
#PAR LAI mean
#p-value: 0.1319
#PAR_LAI_mean_10m<-lm(PAR_LAI_mean~ DEM+tpi10m+lulc_reduced,data=bighornData)
PARLAI.DEM<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=DEM, y=PAR_LAI_mean, color=PAR_LAI_mean)) +
scale_color_viridis_c(option="viridis")+
labs(x = "Elevation", y ="PAR Leaf Area Index Mean",
title ="PAR Leaf Area Index Mean vs Elevation", color ="PAR LAI Mean") +
theme(plot.title = element_text(hjust = 0.5))
PARLAI.DEM
PARLAI.TPI10<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=tpi10m, y=PAR_LAI_mean, color=PAR_LAI_mean)) +
scale_color_viridis_c(option="viridis")+
labs(x = "10m TPI", y ="PAR Leaf Area Index Mean",
title ="PAR Leaf Area Index Mean vs 10m TPI", color ="PAR LAI Mean") +
theme(plot.title = element_text(hjust = 0.5))
PARLAI.TPI10
PARLAI.LULC<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=PAR_LAI_mean, color=PAR_LAI_mean)) +
scale_color_viridis_c(option="viridis")+
labs(x = "Landcover Type", y ="PAR Leaf Area Index Mean",
title ="PAR Leaf Area Index Mean vs Landcover Type", color ="PAR LAI Mean") +
theme(plot.title = element_text(hjust = 0.5))
PARLAI.LULC
PAR_LAI_CAR<-car::avPlots(model = mvPAR_LAI_mean_10m,
col=carPalette()[7],
col.lines=carPalette()[2],
grid=TRUE,
ellipse=FALSE,
marginal.scale=TRUE,
main=("PAR LAI Added-Variable Plots")) 
#PAR LAI sd
#p-value: 0.02881
#PAR_LAI_sd_10m<-lm(PAR_LAI_sd~ eastwest, data=bighornData)
PARLAIsd.EW<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=eastwest, y=PAR_LAI_sd, color=PAR_LAI_sd)) +
scale_color_viridis_c(option="viridis")+
labs(x = "East-West Aspect", y ="PAR Leaf Area Index Sd",
title ="PAR Leaf Area Index Sd vs East-West Aspect", color ="PAR LAI Sd") +
theme(plot.title = element_text(hjust = 0.5))
PARLAIsd.EW
PAR_LAI.sd_CAR<-car::avPlots(model = PAR_LAI_sd_10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("PAR LAI Added-Variable Plot")) 
PAR Average Plots
#PAR Average mean
#p-value: 0.008254
#PAR_Average_mean_10m<-lm(PAR_Average_mean~ lulc_reduced,data=bighornData)
PARAvg.LULC<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=PAR_Average_mean, color=PAR_Average_mean)) +
scale_color_viridis_c(option="magma")+
labs(x = "Landcover Type", y ="PAR Average Mean",
title ="PAR Average Mean vs Landcover Type", color ="PAR Average Mean") +
theme(plot.title = element_text(hjust = 0.5))
PARAvg.LULC
PAR_Avg_CAR<-car::avPlots(model = PAR_Average_mean_10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("PAR Average mean Added-Variable Plots")) 
#PAR Average SD
#p-value: 0.01623
#PAR_Average_sd10m<-lm(PAR_Average_sd~ TRI+slope+ northsouth,data=bighornData)
PARAvg.sd.TRI<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=TRI, y=PAR_Average_sd, color=PAR_Average_sd)) +
scale_color_viridis_c(option="magma")+
labs(x = "Topographic Roughness Index", y ="PAR Average Sd",
title ="PAR Average Sd vs Topographic Roughness Index", color ="PAR Average Sd") +
theme(plot.title = element_text(hjust = 0.5))
PARAvg.sd.TRI
PARAvg.sd.slope<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=slope, y=PAR_Average_sd, color=PAR_Average_sd)) +
scale_color_viridis_c(option="magma")+
labs(x = "Slope Anlge", y ="PAR Average Sd",
title ="PAR Average Sd vs Slope Angle", color ="PAR Average Sd") +
theme(plot.title = element_text(hjust = 0.5))
PARAvg.sd.slope
PARAvg.sd.NS<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=northsouth, y=PAR_Average_sd, color=PAR_Average_sd)) +
scale_color_viridis_c(option="magma")+
labs(x = "North-South Aspect", y ="PAR Average Sd",
title ="PAR Average Sd vs North-South Aspect", color ="PAR Average Sd") +
theme(plot.title = element_text(hjust = 0.5))
PARAvg.sd.NS
PARAvg_Car<-car::avPlots(model = mvPAR_Average_sd10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("PAR Average sd Added-Variable Plots")) 
Canopy Density Plots
#Canopy Density Mean
#p-value: 8.287e-06
#Canopy_Density_mean_10m<-lm(Canopy_Density_mean~ DEM+ lulc_reduced,data=bighornData)
CD.DEM <-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=DEM, y=Canopy_Density_mean, color=Canopy_Density_mean)) +
scale_color_viridis_c(option="rocket")+
labs(x = "Elevation", y ="Canopy Density Mean",
title ="Canopy Density Mean vs Elevation", color ="Canopy Density Mean") +
theme(plot.title = element_text(hjust = 0.5))
CD.DEM 
CD.LULC<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=Canopy_Density_mean, color=Canopy_Density_mean)) +
scale_color_viridis_c(option="rocket")+
labs(x = "Landcover Types", y ="Canopy Density Mean",
title ="Canopy Density Mean vs Landcover Types", color ="Canopy Density Mean") +
theme(plot.title = element_text(hjust = 0.5))
CD.LULC
CD_CAR<-car::avPlots(model = Canopy_Density_mean_10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Canopy Density mean Added-Variable Plots")) 
#Canopy Density SD
#5.175e-05
#Canopy_Density_sd_10m<-lm(Canopy_Density_sd~ DEM+tpi100m+ lulc_reduced, data=bighornData)
CDs.DEM <-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=DEM, y=Canopy_Density_sd, color=Canopy_Density_sd)) +
scale_color_viridis_c(option="rocket")+
labs(x = "Elevation", y ="Canopy Density Sd",
title ="Canopy Density Sd vs Elevation", color ="Canopy Density Sd") +
theme(plot.title = element_text(hjust = 0.5))
CDs.DEM 
CDs.TPI100<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=tpi100m, y=Canopy_Density_sd, color=Canopy_Density_sd)) +
scale_color_viridis_c(option="rocket")+
labs(x = "100m TPI", y ="Canopy Density Sd",
title ="Canopy Density Sd vs 100m TPI", color ="Canopy Density Sd") +
theme(plot.title = element_text(hjust = 0.5))
CDs.TPI100
CDs.LULC<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=Canopy_Density_sd, color=Canopy_Density_sd)) +
scale_color_viridis_c(option="rocket")+
labs(x = "Landcover Types", y ="Canopy Density Sd",
title ="Canopy Density Sd vs Landcover Types", color ="Canopy Density Sd") +
theme(plot.title = element_text(hjust = 0.5))
CDs.LULC
CDs_CAR<-car::avPlots(model = mvCanopy_Density_sd_10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Canopy Density sd Added-Variable Plots")) 
Leaf Angle Plots
#Leaf Angle mean
#p-value: 0.5935
#Leaf_Angle_mean_10<-lm(Leaf_Angle_mean~ lulc_reduced, data=bighornData)
LA.LULC<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=Leaf_Angle_mean, color=Leaf_Angle_mean)) +
scale_color_viridis_c(option="plasma")+
labs(x = "Landcover Types", y ="Leaf Angle Mean",
title ="Leaf Angle Mean vs Landcover Types", color ="Leaf Angle Sd") +
theme(plot.title = element_text(hjust = 0.5))
LA.LULC
LA_CAR<-car::avPlots(model = Leaf_Angle_mean_10,
col=carPalette()[1],
col.lines=carPalette()[1],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Leaf Angle mean Added-Variable Plots")) 
#Leaf Angle sd
#p-value: 0.001394
#Leaf_Angle_sd_10<-lm(Leaf_Angle_sd~ tpi1000m+lulc_reduced,data=bighornData)
LAs.TPI1000<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=tpi1000m, y=Leaf_Angle_sd, color=Leaf_Angle_sd)) +
scale_color_viridis_c(option="plasma")+
labs(x = "1000m TPI", y ="Leaf Angle Sd",
title ="Leaf Angle Sd vs 1000m TPI", color ="Leaf Angle Sd") +
theme(plot.title = element_text(hjust = 0.5))
LAs.TPI1000
LAs.LULC<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=Leaf_Angle_sd, color=Leaf_Angle_sd)) +
scale_color_viridis_c(option="plasma")+
labs(x = "Land Cover Types", y ="Leaf Angle Sd",
title ="Leaf Angle Sd vs Land Cover Types", color ="Leaf Angle Sd") +
theme(plot.title = element_text(hjust = 0.5))
LAs.LULC
LA.sd_CAR<-car::avPlots(model = Leaf_Angle_sd_10,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Leaf Angle sd Added-Variable Plots")) 
Diversity Plots
#Diversity
# p-value: 0.0009821
#Diversity_10m<-lm(div~DEM+northsouth+lulc_reduced, data=bighornData)
div.DEM<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=DEM, y=div, color=div)) +
scale_color_viridis_c(option="turbo")+
labs(x = "Elevation", y ="Diversity",
title ="Diversity vs Elevation", color ="Diversity") +
theme(plot.title = element_text(hjust = 0.5))
div.DEM
div.NS<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=northsouth, y=div, color=div)) +
scale_color_viridis_c(option="turbo")+
labs(x = "North-South Aspect", y ="Diversity",
title ="Diversity vs North-South Aspect", color ="Diversity") +
theme(plot.title = element_text(hjust = 0.5))
div.NS
div.lulc<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=div, color=div)) +
scale_color_viridis_c(option="turbo")+
labs(x = "Landcover Types", y ="Diversity",
title ="Diversity vs Landcover Types", color ="Diversity") +
theme(plot.title = element_text(hjust = 0.5))
div.lulc
div_Car<-car::avPlots(model = Diversity_10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Diversity Added-Variable Plots")) 
DBH Plots
# Average DBH
#p-value: 0.0009843
#DBH_mean_10m<-lm(Average_DBH~ tpi100m+ lulc_reduced,data=bighornData)
DBH_Avg.TPI100<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=tpi100m, y=Average_DBH, color=Average_DBH)) +
scale_color_viridis_c(option="cividis")+
labs(x = "100m TPI", y ="Average DBH",
title ="Average Diameter Breast Height vs 100m TPI", color ="Average DBH") +
theme(plot.title = element_text(hjust = 0.5))
DBH_Avg.TPI100
DBH_Avg.LULC<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=Average_DBH, color=Average_DBH)) +
scale_color_viridis_c(option="cividis")+
labs(x = "Landcover Types", y ="Average DBH",
title ="Average Diameter Breast Height vs Landcover Types", color ="Average DBH") +
theme(plot.title = element_text(hjust = 0.5))
DBH_Avg.LULC
DBH_Avg_CAR<-car::avPlots(model = mvDBH_mean_10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("DBH mean Added-Variable Plots")) 
#SD DBH
#p-value: 0.001216
#DBH_sd_10m<-lm(DBH_Sd~ tpi1000m+ lulc_reduced, data=bighornData)
dbh.sd.tpi1000<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=tpi1000m, y=DBH_Sd, color=DBH_Sd)) +
scale_color_viridis_c(option="cividis")+
labs(x = "1000m TPI", y ="DBH Sd",
title ="Sd of Diameter Breast Height vs 1000m TPI", color ="DBH Sd") +
theme(plot.title = element_text(hjust = 0.5))
dbh.sd.tpi1000
dbh.sd.LULC<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=DBH_Sd, color=DBH_Sd)) +
scale_color_viridis_c(option="cividis")+
labs(x = "Land Cover Types", y ="DBH Sd",
title ="Sd of Diameter Breast Height vs Land Cover Types", color ="DBH Sd") +
theme(plot.title = element_text(hjust = 0.5))
dbh.sd.LULC
lmDBH_sd_10m<-glm(DBH_Sd~ tpi1000m+ lulc_reduced,
family=gaussian(),
data=bighornData)
dbh.sd.CAR<-car::avPlots(model = lmDBH_sd_10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("DBH sd Added-Variable Plots"))
Stand Density Plots
#Stand Density
#p-value: 1.806e-06
#Stand_Density_10m<-lm(Stand_Density~ DEM+ tpi100m+ lulc_reduced, data=bighornData)
den.dem<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=DEM, y=Stand_Density, color=Stand_Density)) +
scale_color_viridis_c(option="mako")+
labs(x = "Elevation", y ="Stand Density",
title ="Stand Density vs Elevation", color ="Stand Density") +
theme(plot.title = element_text(hjust = 0.5))
den.dem
den.tpi100<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=tpi100m, y=Stand_Density, color=Stand_Density)) +
scale_color_viridis_c(option="mako")+
labs(x = "100m TPI", y ="Stand Density",
title ="Average Tree Height vs 100m TPI", color ="Stand Density") +
theme(plot.title = element_text(hjust = 0.5))
den.tpi100
den.lulc<-ggplot(bighornDataNS)+
geom_point(mapping=aes(x=lulc_reduced, y=Stand_Density, color=Stand_Density)) +
scale_color_viridis_c(option="mako")+
labs(x = "Landcover Types", y ="Stand Density",
title ="Stand Density vs Landcover Types", color ="Stand Density") +
theme(plot.title = element_text(hjust = 0.5))
den.lulc
den.car<-car::avPlots(model = Stand_Density_10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Stand Density Added-Variable Plots")) 
Tree Height Plots
#Mean Tree Height
#p-value: 0.0199
#Mean_Tree_Height10m<-lm(Mean_Tree_Height~northsouth+lulc_reduced, data=bighornDataTreeHeightsAndAge)
avgTH.NS<-ggplot(NSbighornDataTreeHeightsAndAge)+
geom_point(mapping=aes(x=northsouth, y=Mean_Tree_Height, color=Mean_Tree_Height)) +
scale_color_viridis_c(option="magma")+
labs(x = "North-South Aspect", y ="Average Tree Height",
title ="Average Tree Height vs North-South Aspect", color ="Average Tree Height") +
theme(plot.title = element_text(hjust = 0.5))
avgTH.NS
avgTH.LULC<-ggplot(NSbighornDataTreeHeightsAndAge)+
geom_point(mapping=aes(x=lulc_reduced, y=Mean_Tree_Height, color=Mean_Tree_Height)) +
scale_color_viridis_c(option="magma")+
labs(x = "Land Cover Types", y ="Average Tree Height",
title ="Average Tree Height vs Land Cover Types", color ="Average Tree Height") +
theme(plot.title = element_text(hjust = 0.5))
avgTH.LULC
avgTH.CAR<-car::avPlots(model = Mean_Tree_Height10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Average Tree Height Added-Variable Plots")) 
#SD Tree Height
#p-value: 0.006275
#SD_Tree_Height10m<-lm(Sd_Tree_Height~ tpi100m+lulc_reduced, data=bighornDataTreeHeightsAndAge)
sd.th.tpi100<-ggplot(NSbighornDataTreeHeightsAndAge)+
geom_point(mapping=aes(x=tpi100m, y=Sd_Tree_Height, color=Sd_Tree_Height)) +
scale_color_viridis_c(option="magma")+
labs(x = "100m TPI", y ="Tree Height Sd",
title ="Tree Height Sd vs 100m TPI", color ="Tree Height Sd") +
theme(plot.title = element_text(hjust = 0.5))
sd.th.tpi100
sd.th.lulc<-ggplot(NSbighornDataTreeHeightsAndAge)+
geom_point(mapping=aes(x=lulc_reduced, y=Sd_Tree_Height, color=Sd_Tree_Height)) +
scale_color_viridis_c(option="magma")+
labs(x = "Land Cover Types", y ="Tree Height Sd",
title ="Tree Height Sd vs Aspect", color ="Tree Height Sd") +
theme(plot.title = element_text(hjust = 0.5))
sd.th.lulc
sd.thCAR<-car::avPlots(model = SD_Tree_Height10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Tree Height sd Added-Variable Plots")) 
#Max Tree Height (northsouth aspect)
#p-value: 0.00969
#Max_Tree_Height10m<-lm(Max_Tree_Height~ aspect+ tpi1000m+lulc_reduced,data=bighornDataTreeHeightsAndAge)
max.th.aspect<-ggplot(NSbighornDataTreeHeightsAndAge)+
geom_point(mapping=aes(x=aspect, y=Max_Tree_Height, color=Max_Tree_Height)) +
scale_color_viridis_c(option="magma")+
labs(x = "Aspect", y ="Max Tree Height",
title ="Maximum Tree Heights vs Aspect", color ="Maximum Tree Height") +
theme(plot.title = element_text(hjust = 0.5))
max.th.aspect
max.th.tpi1000<-ggplot(NSbighornDataTreeHeightsAndAge)+
geom_point(mapping=aes(x=tpi1000m, y=Max_Tree_Height, color=Max_Tree_Height)) +
scale_color_viridis_c(option="magma")+
labs(x = "1000m TPI", y ="Max Tree Height",
title ="Maximum Tree Heights vs 1000m TPI", color ="Maximum Tree Height") +
theme(plot.title = element_text(hjust = 0.5))
max.th.tpi1000
max.th.lulc<-ggplot(NSbighornDataTreeHeightsAndAge)+
geom_point(mapping=aes(x=lulc_reduced, y=Max_Tree_Height, color=Max_Tree_Height)) +
scale_color_viridis_c(option="magma")+
labs(x = "Landcover Types", y ="Max Tree Height",
title ="Maximum Tree Heights vs Landcover Type", color ="Maximum Tree Height") +
theme(plot.title = element_text(hjust = 0.5))
max.th.lulc
max.th.CAR<-car::avPlots(model = Max_Tree_Height10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Maximum Tree Height Added-Variable Plots")) 
Stand Age Plots
#Average Stand Age
#p-value: 0.0001305
#Average_Stand_Age10m<-lm(Average_Stand_Age~aspect+ tpi1000m +lulc_reduced, data=bighornDataTreeHeightsAndAge)
avg.sa.aspect<-ggplot(NSbighornDataTreeHeightsAndAge)+
geom_point(mapping=aes(x=aspect, y=Average_Stand_Age, color=Average_Stand_Age)) +
scale_color_viridis_c(option="turbo")+
labs(x = "aspect", y ="Average_Stand_Age", title ="Average Stand Age vs Aspect", color ="Average_Stand_Age") +
theme(plot.title = element_text(hjust = 0.5))
avg.sa.aspect
avg.sa.tpi1000<-ggplot(NSbighornDataTreeHeightsAndAge)+
geom_point(mapping=aes(x=tpi1000m, y=Average_Stand_Age, color=Average_Stand_Age)) +
scale_color_viridis_c(option="turbo")+
labs(x = "1000m TPI", y ="Average Stand Age", title ="Average Stand Age vs 1000m TPI", color ="Average Stand Age") +
theme(plot.title = element_text(hjust = 0.5))
avg.sa.tpi1000
avg.sa.lulc<-ggplot(NSbighornDataTreeHeightsAndAge)+
geom_point(mapping=aes(x=lulc_reduced, y=Average_Stand_Age, color=Average_Stand_Age)) +
scale_color_viridis_c(option="turbo")+
labs(x = "Landcover Types", y ="Average Stand Age",
title ="Average Stand Age vs Landcover Type", color ="Average Stand Age") +
theme(plot.title = element_text(hjust = 0.5))
avg.sa.lulc
avg.sa.CAR<-car::avPlots(model = Average_Stand_Age10m,
col=carPalette()[3],
col.lines=carPalette()[5],
grid=TRUE,
ellipse=FALSE,
marginal.scale=FALSE,
main=("Estimated Stand Age Added-Variable Plots")) 
ggplot(bighornDataNS)+
geom_point(mapping=aes(x=DEM, y=Average_DBH, color=Stand_Density)) +
scale_color_viridis_c(option="magma")+
labs(x = "Elevation", y ="Average_DBH",
title ="Elevation, DBH and Density", color ="Stand Density") +
theme(plot.title = element_text(hjust = 0.5))